Many-body Landau-Zener tunneling in the Bose-Hubbard model 
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We study a model for ultracold, spinless atoms in quasi-one dimensional optical lattices and 
subjected to a tunable tilting force. Statistical tests are employed to quantitatively characterize the 
spectrum of the Floquet-Bloch operator of the system along the transition from the regular to the 
quantum chaotic regime. Moreover, we perturbatively include the coupling of the one-band model 
to the second energy band. This allows us to study the Landau-Zener interband tunneling within a 
truly many-body description of ultracold atoms. The distributions of the computed tunneling rates 
provide an independent and experimentally accessible signature of the regular-chaotic transition. 
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I. INTRODUCTION 

Bose-Einstein condensates loaded into optical lattices, 
which perfectly realize spatially periodic potentials, rep- 
resent an exciting field of research in the sense that many 
simplified toy models of condensed matter physics can 
be studied in an exceptionally clean manner [T]. This is 
achieved by modern means of atom and quantum optics 
that allow the experimentalist an unprecedented control 
of initial conditions in coordinate and momentum space 
and also of the desired dynamics [U H] . 

A paradigm of quantum transport on the microscopic 
scale is the Wannier-Stark (WS) problem, where parti- 
cles move in a tilted, but otherwise spatially periodic 
potential. The famous Bloch oscillations and related 
phenomena, such as interband tunneling, were observed 
in experiments with quasi-particles in super lattices [3], 
with light in optical nonlinear media [4] , and in great de- 
tail also with ultracold atoms moving in optical lattices 
[H El m El IS] • All experimental studies based on the latter 
realization were performed in a regime where atom-atom 
interactions are either negligible [5] or they reduce to a 
perturbative mean-field effect [51 [TJ |S] . 

The regime of strong correlations in the WS system, 
in which interactions cannot be reduced to a mean-field 
model or even dominate the evolution has been addressed 
only theoretically up to now [3] [ID] HU [12 . State-of-the- 
art experiments are, however, capable of getting into a 
regime of filling factors (i.e. of atoms per lattice site) of 
the order one, where interaction-induced correlations are 
crucial [HIH]. 

Motivated by the experimental progress, we extend 
previous studies of the asymmetric triple-well [TS] and 
of the WS problem UHl [H] ■ In the present work we 
give a more comprehensive and quantitative account of 
our findings briefly reported in |16) . In Section [ll] we in- 
troduce our two-bands Bose-Hubbard (BH) model and 



focus on the dynamics within the first band of the op- 
tical lattice. In contrast to the vast literature which fo- 
cuses on regimes around Mott-Insulator like phase tran- 
sitions in the absence of an additional Stark force, see e.g. 
[T71[TS1 [ini [201 1211 123, we concentrate on the BH model 
in the superfluid realm and in the presence of a static tilt. 
We characterize the transition between the regular and 
the chaotic realm of the quantum spectra by a quantita- 
tive and systematical analysis based on statistical tests. 
In Section [ni| we perturbatively include the decay to the 
second band via Landau-Zener like tunneling processes. 
The quantum spectrum of the latter, non-unitary prob- 
lem is analyzed in Section |IV| and found to essentially 
reproduce the properties of the purely one-band approx- 
imation down to rather small lattice depths. The re- 
sulting decay rates for the interband tunneling strongly 
depend on the many-particle nature of the problem and 
are found to correlate with the transition to the quantum 
chaotic regime. As a consequence, signatures of many- 
body quantum chaos are predicted to be accessible to 
experiments with ultracold atoms over a broad range of 
parameters, in both "horizontal" transport along the lat- 
tice and in interband "vertical" transport. Our results 
are finally summarized in Section [V] 



II. SPECTRAL ANALYSIS OF THE 
ONE-BAND BOSE-HUBBARD MODEL 



We briefiy review the general Hamiltonian for a system 
of spinless, interacting atoms in a quasi-one dimensional 
optical lattice subjected to a tunable tilting force F. We 
start out with the purely periodic problem, = 0, in an 
optical potential of spacing a, recoil momentum — 7r/a 
and typical kinetic energy En — fci/2m. The optical 
potential and the kinetic energy form the single-particle 
Bloch Hamiltonian: 
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The eigenfunctions are the Bloch waves -01 labeled by 
the quasimomentum k and the band index a [23l . with 
dispersion law E^ a- In the Appendix [X| we explain how 
to use the single-particle solutions ipk.a to build a set of 
localized orbitals Xi,a called Wannier Functions (WF). 
hi the limit of deep lattices, the orbital xe,a goes to the 
wave function of the a— th excited level for an harmonic 
potential centered on the ^— lattice site. We use the WF 
to expand the field of the ultracold bosons: 



(2) 



Then we introduce a Stark force F that tilts the optical 
potential and a zero-range interaction between the atoms, 
parametrized by the scattering length as- In a quasi-one 
dimensional optical lattice - as realizable in experiments 
[l4| - the scattering length is derived from the true three- 
dimensional scattering length via a renormalization that 
accounts for the transverse confinement of the atomic 
wave functions |14j and the physical dimension is then 
[as] = L~^. The Hamiltonian in the second quantization 
is written as: 

H = J ^H^) [i?o„o(x) + Fx] 4>{x)dx 

liiras /■0t(a.)0t(a.)0(a.)0(a.)da;. (3) 

Substituting the expansion Eq. ^ into Eq. ^ we obtain 
the Hamiltonian in terms of the creation and annihilation 
operators ^ , dg^a i for a particle in the £—th site and 
the a-th energy band of the lattice. The number oper- 
ator is n£^Q — d\^ag^a- We restrict the analysis to the 
first two bands, which can be addressed by experiments 
[24j and that can be handled numerically without great 
difficulties. The coefficients of the Hamiltonian are given 
by integrals involving the WF: the exact computation of 
the WF outlined in Appendix|X]motivates the selection of 
the operators that are most relevant for V > En. We are 



left with the on-site energy a\ ^ai^ai the kinetic energy 
'^^+1 ofli.ai and the on-site interaction between atoms in 
the same band fig^aifii.a — 1) ^5^, for a e {1,2}. More- 
over we have on-site interaction between atoms in differ- 
ent bands ng ini 2 and two transition operators 0^2^^,!! 



i^f, 20^.2 that are the subject of a detailed analysis 



in Section Hill 

The dimension Z^h of the Hilbert space spanned by the 
Fock states |n) (defined in Appendix [b|) , for a system of 
N bosons distributed over L lattice sites, occupying up 
to the 2nd band of the periodic potential, is given by the 
combinatorial formula Dii = {N + 2L- 1)\/N\{2L ~ 1)!. 
The typical number of lattice sites in experiments is 
L < 100 and the filling-factor N/L is of order unity 
[131 Hj, such that the exponential increase of Z3h with 
the system size limits any exact numerical approach to 
smaller systems, where we impose the cyclic boundary 
conditions iiL.a = cio,a- The implementation of these 
conditions requires the system to be translationally in- 
variant. The Stark potential Fx, however, spoils the 
periodicity of the Bloch Hamiltonian Eq. Q and pro- 
duces localized WS eigenstates instead of traveling Bloch 
waves [23]. We follow [10 and proceed to eliminate the 
Stark potential from the Hamiltonian by changing to the 
Interaction Representation (IR) with respect to Hs — 
F J (f)^x(f>dx = aF J2e=i^a^^^,a- The Hamiltonian in 
the IR, Hit) = e~'"''*He+^"'^\ is time-dependent, and 
the problem becomes conceptually more complicated. 

We rescale the energies by Ej^, the lengths by a, 
the momenta by and we make the substitutions 
F ^ FEjf/a, X ^ x/V^- The on-site energi es Eg 
and the hopping amplitudes are given in Eq. ( A3 1 . 
The interaction coefficients Wa, Wx, are proportional to 
the coupling constant W = Anas /amEn and are given 
in Eq. (|A5[). The "dipole" coefficient dp is given in 



Eq. (|A6|. The Hamiltonian of Eq. restricted to the 
first two bands of the lattice, finally reads in the IR: 



* ^ f 1 1 

H{t) = J2 ] - 2'^i<^'^'"-l+i.i^'iA + H.c. + -Wifie^iifiLi - 1) + (1 ^ 2) 

e=i ^ 

+ 2Wxne^ihi^2 + Fdp{al.^de^i + il.c.) + ^Wx{al^^dl.^di^2d£.2+ H.c.) 



(4) 



In the IR the Hamiltonian is again symmetric for discrete 
translations in space and it has lost the time indepen- 
dence but it is periodic with the Bloch period Tb = 2-k/F. 
We assume that the initial state, for = 0, is super- 
fluid, characterized by the delocalization of the atoms 
over the entire lattice. The critical conditions on the pa- 
rameters, that enforce the superfluid phase in the present 



context [T7] is Wx/ Ji < 5.8. Following [9 , in the present 
Section we set the lattice depth ~ 5, which gives 
Ji ~ 0.038 and the interaction coefficient W ~ 0.016, 
with Wi ~ 0.032. 

The object of the subsequent study is the evolution op- 
erator up to the Bloch period Z^fb = exp(— J^^ iH(t)dt), 
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FIG. 1: The quasienergy levels E of the FB operator for a 
system with A'' = 4, L = 5, in the subspace with k, = Q. The 
WS ladder is seen on the left panel, the perturbative splitting 
of the first rung is magnified in the lower panel. In the upper 
panel a single linear function of 1 /F is added to 2E/F to elim- 
inate an overall winding trend, and allow a better visibility of 
the avoided crossings. The parameters of the Hamiltoniana 
are Ji = 0.038, Wi = 0.032, as used also in the subsequent 
Figures. 



called Floquet-Bloch (FB) operator ^lOJ. The results pre- 
sented in the following confirm and extend the results of 
[lOj . The discrete translational symmetry of the Hamil- 
tonian entails that the FB operator is a block-diagonal 
matrix in the basis of Eq. ( |B2[ ) , labeled by a many-body 
quasimomentum k. The dynamics of the atoms in the 
lattice is complex because many vectors take part in the 
time evolution of an arbitrary initial state. The strong 
mixing of the basis vectors in time means that the evolu- 
tion of a state is not bound to a small subspace of the to- 
tal Hilbert space (contrary to [E]), but, after a Bloch pe- 
riod, the initial state spreads over the entire Hilbert sub- 
space with definite quasimomentum (e.g. k = 0). This is 
evidenced by the dependence on F of the quasienergies 
Ej, obtained from the eigenvalues exp(— iE'jTB) of the 
FB operator. 

In the Fig. [l] we show the quasienergy spectrum as the 
"control parameter" F is varied. In the limit F — > 0, 
we recover the standard BH model (the analysis of the 
FB spectrum is here, however, not useful since the Bloch 
period tends to infinity) . In the regime where the atomic 
interactions are negligible with respect to the lattice po- 
tential, the single particle Bloch picture is adequate and 
the spectrum is simply a finite band. For F > 0.1 the 
single-particle WS ladder is found, i.e. a fan of energy 
levels Em{F) ~ 27rmF, m integer. Since the interactions 
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FIG. 2: The distribution P(s) and the CDF (insets) for the 
quasienergy spacings (stairs), the WD (solid) and the Poisson 
distribution (dashed). The parameters are A'' = 5, L = 8, 
logj^g(27r/F) ~ 2.0 (left, regular) and 2.4 (right, chaotic). 



are non-zero, the levels are split up and the first order 
perturbative effect on the ladder was computed in |llj . 
The central WS rung is split into levels which are pro- 
portional to the interaction energies of many atoms in 
a site, Win{n — 1), n = 1, . . . ,4. Since the level split- 
tings have the common factor Wi , a collapse and revival 
of quasimomentum oscillations (the Fock space version 
of the single-particle Bloch Oscillations) was predicted 
for this parameter range in [11 . On the contrary, an ir- 
reversible decay of the quasimomentum oscillations was 
found in the range F ~ Wi 9 , where avoided cross- 
ings dominate the spectrum. These are directly, exper- 
imentally observable consequences of the complex level 
structure presented in Fig. [T] The presence of avoided 
crossings means that a strong mixing of the Fock states 
is necessary to build the eigenstates of the system, and 
no set of quantum numbers can be assigned to individ- 
ual levels as F varies. In the region of parameters where 
the energy scales of the system are comparable in mag- 
nitude, we can characterize the spectrum in terms of the 
statistical distribution of the quasienergies spacings and 
by statistical measures for the eigenfunctions. The latter 
have been intensively studied in |26j . In the following, we 
concentrate on the statistical behaviour of the quasiener- 
gies, which is closely linked to the behaviour of the open 
system studied in Section |III[ 

The probability P(s)ds that the magnitude of a given 
interval spacing Sj — AEj/{AEj)j is in [s, s+ds] is given 
by the Poisson distribution P(s) = exp(— s) |27j for an 
uncorrelated spectrum in the regular regime. Strongly 
correlated quantum spectra, corresponding to the chaotic 
regime in our many-particle model, are well modeled 
by the Wigner-Dyson (WD) distribution for a Circular 
Orhogonal Ensemble of random matrices [27/. P(s) = 
TT s exp(— TT /4)/2. In Fig. [2] the probability distribu- 
tion and the Cumulative Distribution Function (CDF) of 
the quasienergy spacings are shown for two paradigmatic 
values of F. The presence of avoided crossings in the 
chaotic case log]^Q(27r/F) ~ 2.4 shows up as a depletion 
of small quasienergy spacings, and the probability to find 
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a level crossing vanishes. 

We improved the statistical description of the 
quasienergy spectrum with further analyses, shown in 
Fig. |3] In the panel (a) we quantify the convergence of 
the quasienergy spacings distribution to the WD profile, 
thus filling the gap between the two pictures of Fig. [2] 
We computed the FB operator for several values of F 
and confronted each spectrum with the WD distribution 
using a modified test, computed as follows. Each 
sequence of levels spacings was algorithmically binned 
to leave 5 < Ob < 10 "observed" spacings in each bin 
b— 1, . . . , iVf, |55]. The "expected" values Ef, are the in- 
tegrals of the theoretical distributions over the bins and 
the sum Q — J2b i^b — Eb)'^ / Et, was calculated. The 
values of Q are distributed according to a distribu- 
tion with iVf, — 1 degrees of freedom and mean N(, — I. 
The renormalized variable 

X^=logio[g/(iV,-l))] (5) 

is thus appropriate to compare several data sets, each 
binned optimally and independently. For F < 0.025, 
is in the range [—0.5, 0.5] (in the bulk of the original 
distributions before the transformation of Eq. ([5| was 
performed), and the correspondent distribution of the 
spacings is well described by a WD profile. As the exter- 
nal force increases and the parameter 2t:/F diminishes, 
the larger values of x^ (lying in the tails of the original x^ 
distributions) indicate that the spectrum is not well char- 
acterized by a WD distribution. The condition for the 
regular-chaotic transition can be directly read off from 
the quantitative statistical test results and corresponds, 
e.g. to logio(2^/f^) ~ 2.3. 

We found that the profile of the quasienergy spacings 
distribution changes smoothly, but we used two statisti- 
cal tests introduced in [2S] (eq. (27) and (29) therein) to 
show that the Poisson and WD statistics are clearly iden- 
tified at the borders of the transition. Fig.|3](b) shows the 
explicit results for the T function of [29 whose predicted 
linear scaling T ~ In s for the Poisson and T ^ 2 In s 
for the WD expectation is confirmed. In Fig. |3] (c) we 
show the variance S^(d£') \W of the number of levels 
N{E, E + dE) found in a finite energy interval dE: 

T.'^idE) = {[N{E,E + dE) - NdEf)E, (6) 

where the average is taken over all the energies and we 
rescaled the spectrum so that the average number of lev- 
els per unit of energy N is equal to unity. A linear and 
logarithmic scaling is predicted for a Poisson and WD- 
like spectra, respectively. The logarithmic behaviour 
clearly prevails over all energy ranges for the chaotic spec- 
trum, only apart from oscillations which are indeed typ- 
ical for samples of finite size (see ^30j for details). 
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FIG. 3: Statistical tests applied to the system of Fig. [2] (a) 
X^— like test of Eq. ([5|. (b) T test and (c) variance of number 
of levels (with mean spacing normalized to 1), for the chaotic 
(circles) and the regular (pyramids) case. In both panels the 
solid and dot-dashed lines are the theoretical predictions for 
the WD and Poisson statistics of energy levels, respectively. 



III. INTERBAND COUPLING IN THE 
MANY-BODY REGIME 

We develop a perturbative analysis of the two-bands 
system of Eq. Q . We consider the many-body dynamics 
within the ground band and the perturbative action of 
the operators that couple the Fock states of the ground 
band to states in the second band. As a consequence of 
the interband terms each Fock state |n) (see Appendix [b| 
suffers an energy shift 6E(n) — iT-p(n)/2, where Tp(n) is 
its decay width. The Hamiltonian is modified accordingly 
and becomes a non-Hermitian effective Hamiltonian for 
the ground band, that yields a non unitary FB operator. 
In the following, we compute the set of decay widths. In 
the next Section we then study the spectrum of this new 
FB operator. 

We first define a basis of unperturbed states. We 
choose to neglect the hopping in the lower band, where 
the WF are more strongly localized, so an unperturbed 
state projected on the Hilbert space of N particles in 
the ground band is a Fock state \n;N). In the second 
band we neglect the interactions, since in the perturba- 
tive approach only a few particles (one or two in the fol- 
lowing) populate the excited levels. So an unperturbed 
state projected on the second band is the solution of the 
one-particle WS problem [THl [5^, i.e. a localized wave 
function centred at site w, written with the Bessel func- 
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tion Jm{x) as: 



+00 



vac) 



(7) 



We approximate the Hilbert space of the system as the 
tensor product of the spaces of the two bands and the 
entanglement between the ground and the excited par- 
ticles is neglected. Then an unperturbed state with one 
or two promoted particles is of the form \n; N — 1) ^ jw) 
or \n; iV — 2) (g) \w, w'), respectively. In the following, the 
occupation number n£^i for the i—th site in the ground 
band is written n^. 

The Hamiltonian Eq. Q contains two mechanisms 
that promote particles to the second band. The first is a 
single-particle effect, a consequence of the external force, 
proportional to the dipole coupling dp between the WF 
of different bands. The Hamiltonian of the perturbation 
is 



Hi = FdF^(a 



(8) 



The second perturbation is a many-body effect, describ- 
ing two particles of the first band that collide and trans- 
form their interaction energy into kinetic energy, entering 
the second band together: 



1 



(9) 



The expectation value of Hi, H2 on the unperturbed 
states, equal to the first-order energy shift 5E(n), is zero 
because the operators do not conserve the number of par- 
ticles Ua within the bands. 

Let us focus on Hi and compute its matrix element for 
the channel: 

\n;N) (g) |vac) -> \n'\N -\) ® \w), < = nu - 1. (10) 

The decay width at first-order is given by the Fermi's 
Golden Rule and only the first term in Eq. (|8| gives 



nonzero contribution for the channel of Eq. (10 1. Our 
result for the matrix element is: 



(/e|(n'|^(a|2a£,i)|n)|vac) 



(11) 



j{-J2/F) S{n'^,ni- I) y/iM J]^ (5(7iJ„,n„). 



The Kronecker 6{-, •) functions act as a selection rule for 
the Fock states that are coupled by the perturbation. 
The tunneling mechanism does not include any income 
of energy from an external source, so the initial and final 
energies, 

Eo{n) = (vac|(n|_ffo|"^)|vac), 
Eoin',w) = {w\{n'\Ho\n')\w), (12) 



must be equal as required by the Golden Rule. The con- 
dition on the energy conservation is relaxed to account 
for the uncertainty AE{n) of the unperturbed energy lev- 
els of the initial and final states. The energy uncertainty 
and the level density function p{E, n) are derived from 
the perturbative action of the hopping operator of the 
first band that has been neglected so far. We postpone 
the computation of these quantities to the end of the 
present Section. The relaxed energy conservation rule se- 
lects from Eq. (Ill the set K of permitted decay channels 



(/i, w) parametrized by the two indices h, w such that: 



Eo{n\w)- Eo{n) - 

= £2 - £1 - F{h ~W)~ Wl{72h " 1) 

AE{n) + AE{n') AE{n) + AE{n') 



.(13) 



The last equation means that the energy S2 — £i required 
to promote a particle to the second band is supplied by 
the decrease of the repulsion energy (proportional to W) 
and by the work of the force (proportional to F) exerted 
on the promoted particle. 

The total width Ti{n) for the decay via the allowed 
channels K, is proportional to the square of the matrix 
element and to the level density p{E, n) given below in 
Eq. (p2|). We arrive at 



ri(n) = 277 F^(fp 

X {\Jh-U-J2/F)-V^i 

{h.w) e K ^ 



(14) 



1 



AE{n)AE{n') 



The perturbation H2 is treated in a similar way, with the 
difference that two particles are promoted to the second 
band, and the position of the second Stark state is 
an additional degree of freedom for the transition. The 
decay channels are: 



\n,N) (g) |vac) |r7;',iV - 2) (g) \w,w'), 



nh 



(15) 



The approximate energy matching equation selects a set 
K of possible decay channels, parameterized by the three 
site indices {h,w,w'): 



{h,w,w') G K s.t. Et){fi' ,w,w') — Et){fi) = 
= 2{s2 - £1) - F{2h -w-w')~ Wi{2nh - 3) 
AE{n) + AE{n') AE{n) + AE{n') 



2 ' 2 

We state the result for the decay width: 

T2{n) = Kw^ Y l\Jk-U~J2/F) 

(h.w.w') e K 

X Jh^^,{-J2/F)\^ Afihirih-l) 



(16) 



AE{n)AE{n') 



(17) 



With respect to Eq. (14 1, the additional degree of free- 



dom w' results in an extra summation extended over the 
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(infinite) possible values of the difference w — w'. This 
follows from the possibility to conserve the energy even 
if a particle is pushed very far, if the other particle is 
pushed almost equally far in the opposite direction. Since 
the decay width Eq. ([T7| depends on the product of two 
(rapidly decaying) Bessel functions we apply the trunca- 
tion — w'l < \J2/F\, to reduce the formula to a finite 
form. 

Now we conclude the computation and derive the en- 
ergy broadening AE{n) of the Fock states in the ground 
band necessary to implement the conditions of Eq. (13) 
and (16 1. In the ground band, the unperturbed Hamil- 



tonian consists only of the on-site interaction operator: 



Ho 



In the case of a single particle in a periodic potential 
the use of first order perturbation theory is wrong, as 
the second order of the perturbation theory diverges be- 
cause of the exact degeneracy in energy of neighbouring 
sites, entailed by the translational symmetry of the lat- 
tice. On the contrary, in the present case, the unper- 
turbed Hamiltonian is just a rough approximation of the 
true Hamiltonian of the system and the remaining oper- 
ators are supposed to remove the degeneracies, since the 
translational symmetry is broken by the external field in 
the WS picture. The perturbation Hamiltonian is given 

by 



H, 



(18) 



and its matrix elements between Fock states are 
1 



{n'\Hh\n) = 



L 



e=i Ae=±i 



n 



xS{n'f^ -i^,ne^i - 1) (5(n^+A£4> "•^+A;a + !)■ (19) 

Transitions are allowed between Fock states that differ 
for one boson in two adjacent holes m, m + Am. The 
transition channels are written as |n) 

nm,l - 1 and n^j+Am.l = "m+Am,l 

the condition of the energy conservation 



\n'), with n'^ i = 
1, and must fulfil 



Eo{n') - Ea{n) - Wi {n^+AmS - n„,,i + 1) = 0. (20) 

The total uncertainty is obtained by adding up the con- 
tributions from all the transitions and the summation 
over the allowed channels can be re-casted in a summa- 
tion over the lattice sites. Using the Golden Rule, each 
Fock state is attributed the following energy uncertainty: 

AE{n) = ^nJl ^ AE{n n') = 

rt' 

= ^ttJi ^ ^ <5(n<>+A^i + l,n£4). (21) 

I Al=±l 




FIG. 4: (a,b) the probability distribution P and (c,d) the 
CDF for the logarithm of the decay widths V in the regu- 
lar regime [(a,c) with logj^Q(27r/_F') ~ 1.5] and in the chaotic 
regime [(b,d) with logj^Q(27r/_F) ~ 2.1]. The size of the system 
IS N = S, L = 7 . The dashed line is the fit with a lognormal 
distribution. The inset in panel (d) shows that the lognormal 
is not appropriate to fit the tails of the distributions in the 
chaotic regime. Here and in the following Figures, V — 1.5, 
£2 - £i = 2.63, Ji = 0.22, Ja = -1.0, Wi = 0.2, = 0.1, 
dp = -0.2. 



Finally, the level density p{E, n) around the unperturbed 
energy Eo{n) of a Fock states \n) is approximated by a 
rectangular profile, of width AE(n) and area unity: 

piE,n) = x{\E - Eoin)\ < AEin)/2} / AE{n). (22) 



IV. RESULTS OF THE PERTURBATIVE 
OPENING OF THE ONE-BAND SYSTEM 

The total width Fp = Fi -I- F2 is now added to the 
single-band Hamiltonian as a complex shift. Given the 
translational symmetry of this Hamiltonian, we added 
the shift to the diagonal in the representation of the cyclic 
basis \an) in Eq. (B2|, as — iFF(cr)/2. The eigenvalues of 
the FB operator are no longer unitary and the quasiener- 
gies Ej has a complex part —iTj/2. We analyzed the de- 
cay rates F^ along with the quasienergy spacings statis- 
tics to study how the dynamics within the first band 
infiuences the coupling to the second band. We first re- 
ported the distribution P(F) in [16j . Here we refine our 



analysis and we first focus on P(A), with A = log^gF, 
shown in Fig.|4]for some paradigmatic cases. The widths 
are much smaller than unity, consistently with a pertur- 
bative approach of the system, yet the lattice potential 
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FIG. 5: (Upper panel) the x —like test for the quasienergy 
spacings s and (lower panel) the spread of the distribution 
of the decay widths F. The size of the system is A'' = 8, 
L = 7. The dashed lines are a guide to the eye and suggest 
that the transition to the chaotic regime can be appreciated 
by looking at both, the real and the imaginary parts of the 
FB eigenvalues. 



is only V = l.5Ej^ to increase the spread of the Stark 
state in the second band. Moreover, the decay channels 
are activated by an increase of the interaction energy, 
which can be experimentally achieved by acting on the 
transversal confining potential |31j of the quasi-one di- 
mensional lattice, or by a Feshbach resonance [T3]. In 
this Section, W ~ 0.02 used in |9j is multiplied by a 
factor of order 10, a value that is still well within the 
experimental possibilities [m |32] . 

In Fig. |4] we compare the distribution of the decay 
widths for two values of F that belong to the reg- 
ular (a,c) and the chaotic region (b,d). The differ- 
ence in the average decay width (A) is due to the im- 
proved energy matching provided by a stronger exter- 
nal field F, that supplies the necessary energy to pro- 
mote particles to the second band. For the parameters of 
Fig.ra the single-particle Landau-Zener formula [23] gives 
Flz = F/(27r)exp[-^2(^^ _ e,)y{8F)] ~ lO'^o, lO'^^ 

for (a,b). Then we see that there are regimes where the 
many-body interactions affect substantially the single- 
particle tunneling rates and cannot be neglected. More- 
over, even mean-field treatments of the Landau-Zener 
tunneling [HI I33j cannot account for several decay chan- 
nels. 

The bulks of the distributions, both in the regular and 
in the chaotic regime, are appropriately fit by a lognormal 
profile P(A) « exp{-[ln(A-A^i„)]V2AA2}/(A-A^i„). 



FIG. 6: (Color online) (a-c) l.h.s of Eq. (231 in the regular 



(a) and in the chaotic regime (b,c), for a system with N — 8, 
L — 7. The dashed area in panels (b,c) shows the part of the 
histogram where the total amount of points is about 40, less 
then 10% of the full sample . The red solid line is the linear 
fit to the profile in the scaling region, (d) The exponent of 
the power law P(r) ~ F~" as obtained from the linear fits 
shown in (a-c), for N — 8, L = 7 (black circles), N — 8, 
L = 5 (red squares), TV = 7, L = 8 (blue diamonds), N = 7, 
L — 5 (green pyramids). To show together different data sets 
we horizontally translated the first by +0.1 and the last by 
-0.05. 



We notice that the spread AA is reduced in the chaotic 
case, where the Fock states are strongly mixed by the 
dynamics and a fast decaying Fock state (in the \aK) 
representation) could be a privileged decay channel for 
many eigenstates of the system. Many eigenstates then 
shared similar decay widths and their statistical distri- 
bution would be thinner. Following this reasoning, we 
interpret the thinner distributions found in the panels 
(b,d) of Fig. |4] as a signature of the strongly correlated 
dynamics. 

This picture is supported by Fig. [5] where the de- 
pendence on F of the spread AA (lower panel) is con- 
fronted with the regular-chaotic transition evidenced by 
the x^— like test of Eq. ^ (upper panel). The shrinking 
of the decay widths distribution goes along with the tran- 
sition, though the precise determination of a transition 
point is precluded by a substantial amount of noise. The 
finite size of the samples that can be managed numeri- 
cally accounts for the noise, as we verified that the profile 
becomes sharper while increasing the size of the system. 
Moreover, since the average decay width decreases with 
smaller F, we need a more precise and hence a more time- 
expensive numerical computation to determine AA as we 
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enter the chaotic regime. 

Finally, we found in [l6j that the tails of the distribu- 
tions in the chaotic case follow the expected power-law for 
the diffusive regime of an open quantum chaotic system. 
In our case, the opening of the ground band subsystem 
is defined by the interband coupling, which in a sense 
attaches "leads" to aZZlattice sites within the sample. In- 
deed, the inset of Fig. |4 (d) shows that the lognormal is 
not a good fit to the tails in the chaotic regime: the dis- 
tributions transform to a power-law P(r) ~ F"" in close 
analogy to the transition from Anderson-localized to dif- 
fusive dynamics in open disordered systems [311 [35] . In 
particular we found a ~ 2, in accordance with the general 
results of [36]. In the Fig. [6] panel (d) we show that this 
value is indeed peculiar to the chaotic regime, and strong 
fluctuations with F mean that the exponent a is not de- 
flned within the transition region. Here we evaluate the 
integrated profile 



1 - CDF(r) - 1 - / P(r')dr' - r 



(23) 



and fit the latter to reduce the statistical fiuctuations due 
to the finiteness of our samples. The finite-size effects are 
more relevant in the chaotic case (panels b,c), because the 
reduced spread of the distribution makes the fit sensible 
to the few points that fall in the further part of the tail 
(shaded in the pictures) . The largest analyzed sample has 
only ~ 400 energy levels, yet the uncertainty on a ~ 2 
on the chaotic side of the transition is less than 10%. 



V. CONCLUSIONS 

We studied the problem of a many-body atomic system 
in a tilted optical lattice, using a statistical analysis of 
the complete quantum spectrum. We extended previous 
work [9j [lOl [16] and verified thoroughly the transition 
from a regular to a quantum chaotic spectrum, found 
in [lOj as the Stark field is varied. The transition takes 
place when the tunneling amplitude in the lattice be- 
comes comparable to the interaction energy of the atoms. 
In this regime it is not possible to find a set of quantum 
numbers that decompose the spectrum into subspaces not 
mixed by a change of the Stark field. Because of the 
strong mixing of the energy eigenstates, many avoided 
crossings are found in the spectrum and the statistical 
analysis of the energy spacings show a characteristic de- 
pletion of small spacings - a signature of quantum chaos. 

We derived a reasonable two-bands model, on the ba- 
sis of which we opened the ground band subsystem by a 
perturbative coupling of the Fock states to excited Stark 
states. We obtained the set of decay rates from the 
complex- valued quasienergies of the FB operator. We 
analyzed the real part of the quasienergies and verified 
that the transition from the regular to the chaotic regime 
is still visible and not much modified with respect to the 
one-band system. Moreover we analyzed the statistical 
distribution of the imaginary part of the quasienergies. 



i.e. the decay rates of the states in the ground band. We 
found that the distributions of the decay rates become 
thinner as the regular-chaotic transition is crossed. We 
reckon that thinner distributions of the decay rates are a 
signature of the complex dynamics, where a few strongly 
decaying states act as leading decay channels. The statis- 
tical characterization of the tunneling rates could be used 
to compute the expected atomic current from the ground 
band to the excited band of the lattice, thus providing 
an experimental probe for the regular-chaotic transition. 
Time-dependent observables could possibly be computed 
also with advanced mean-field techniques as reported in 
Ref. [37], as long as no full spectral information of the 
many-body system is wanted. 

Of course, a deeper investigation is desirable to under- 
stand on quantitative grounds the full interplay of the 
dynamics within the first band and the decay towards 
higher bands, and a more detailed analysis of the inter- 
band coupling will be worthwhile in a full-blown model in 
which at least two bands are completely included. Our re- 
sults are a first step in this direction of studies of regimes 
in which both "horizontal" and "vertical" quantum trans- 
port are simultaneously present and influence each other 
in a complex manner. 
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APPENDIX A: COMPUTATION OF THE 
COEFFICIENTS OF THE MANY-BODY 
HAMILTONIAN 

The coefficients of the Hamiltonian Eq. Q are given 
by integrals of the WF xe a , once the expansion of Eq. ^ 
is substituted into Eq. ([3|. The WF are defined as: 



2tt ^ 

e 



(Al) 



From the definition it follows that Xi.ai^) — Xai^ — £a), 
so that all the Bloch waves for a band can be computed 
starting from a single Wannier function Xa- We com- 
puted the WF following the method introduced in |38j . 
A Bloch wave is factorized as ipk,a{x) = e^''^Uk,a, where 
the periodic function Uk,a{x) = Uk,a{x + i) is expanded 



9 



in a truncated basis of momenta using multiples of 27rj/a 
only: 



Q-i 



(A2) 



"0,0 = n'L,a allow us to define a shift operator S [TU] 
that translates the occupation numbers of a Fock state 



,nL- 



(Bl) 



Using a gauge transform (p — > — i^a; + an effective 
Hamiltonian for the periodic Uk,a is derived from Eq. ([ij. 
We solved the effective Schrodinger equation as a linear 
system for Uj{k,a) with dispersion law Ek^a as eigenval- 
ues. We obtained the on-site energies and the hopping 
amplitudes from the Fourier transform of Ek^a'- 



+1 



Ek,a'ik, Ja — 



r + 1 

J Ek^^e'^'^dk. (A3) 



The WF is finally computed from the inversion of 
Eq. (Al I. The relative phase of the eigenvectors Uj{k, a), 



for different quasimomenta, is not unique. For the simple 
case of a sinusoidal lattice potential, the correct choice 
of the phases can be inferred from [55] : 

u,{k,l) ^ \uj{k,l)\, uj{k,2) |uj(fc,2)|sign(2j + fc). 

(A4) 

This phase choice guarantees the correct inversion sym- 
metry Xa{—x) = {—l)°'~'^Xa{x) of the WF of the first 
and of the second band, respectively. The interaction 
coefficients for Eq. Q read: 



Wa = W 



/OO nOQ 
Xl dx, W^=W xlxl dx. (A5) 
-OO J —OO 



The WF are an orthogonal set of functions, so that differ- 
ent bands are decoupled in the one-body dynamics. The 
tilting potential Fx has nonvanishing "dipole" matrix el- 
ements that couples adjacent bands: 



Xi{x)xx2{x) dx. (A6) 



APPENDIX B: COMPUTATION OF THE 
FLOQUET-BLOCH OPERATOR 

The Fock states \n) are defined by the sequence of oc- 
cupation numbers {ng^f^-^ for the L sites of the lat- 
tice, in each band a. The cyclic boundary conditions. 



The operator S naturally decomposes the Fock space into 
equivalence classes of vectors generated by repeated ap- 
plication of S onto a "seed" vector |cr) with M{a) < L 
such that = \n) for all the vectors \n) in the 

class. A new basis Ictk) can be introduced for which 
the many-body quasimomentum k — j /M(a) (0 < j < 
M((t)), supplies a convenient label: 



M{a) 



(B2) 



The cyclic basis decomposes the Hamiltonian Eq. Q into 
a block form that transfers to the FB operator, whence 
Ufb = ffij^i ^fb(k = i/L) with the obvious advantage 
that we can diagonalize separately each block of dimen- 
sion D < Dn/L. This decomposition not only leads to 
a substantial numerical simplification but is also of dy- 
namical relevance [TIJl. Moreover, we exploited that the 
Hamiltonian matrix is sparse and we verified that the 
fraction of the nonzero entries is 4.0 x D^^-^ in the limit 
of large Hilbert spaces D ^ 1. 

The column c of the FB operator coincides with the 
column of the coefficients of the basis state of index 
c, evolved in time up to one Bloch period. We used 
a fourth-order Runge-Kutta time integrator with adap- 
tive stepsize, tuned in precision by the upper bound e 
of the estimated one-step error [2F. The value of e 
was chosen to suppress, up to Tb, the well-known ex- 
ponential instability of the Runge-Kutta method applied 
to the Schrodinger equation. The quantity Q{E, E') = 

1/2 

(Si ^ -^iP/-^^) was used to compare different 
spectra {E}, {E'} and we verified the consistency of 
our computations, finding a power-law self-convergence 
g(i;w,£;('-i)) cx 4-2 for a sequence of tests with in- 
creasing required precision, {e^jr 0. The achieved 
precision scales with cpu time t as Q cx t~^'^ . A preci- 
sion up to 10^^^ was necessary to reliably compute the 
small complex part in the eigenvalues of the FB operator. 
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